Bioconductor

title: 'RGMQL Example R Notebook: Use case 3' author: "Silvia Cascianelli" date: "r Sys.Date()" output: pdf_document: default html_notebook: default html_document: default BiocStyle::html_document: chunk_output_type: inline

#knitr::opts_chunk$set(cache = TRUE)

In this example we investigate the ENCODE Chip-Seq narrow and broad data, extracted from the GMQL remote repository to automate all the steps needed to identify transcription factor (TF) high accumulation DNA zones using RGQML together with TFHAZ, another R/Bioconductor package. The knowledge of DNA regions in which transcription factors bind, in particular the HOT (High Occupancy Target) regions occupied by many different factors, is crucial to understand cancer genesis and develop new targeted therapies.

Load the RGMQL package and initialize the remote GMQL context of scalable data management engine, specifying remote_processing = TRUE, and, possibly, an authenticated login:

library(RGMQL)
remote_url <- "http://www.gmql.eu/gmql-rest"
init_gmql(url = remote_url, remote_processing = TRUE) #, username = 'XXXX', password = 'XXXX')

Download and extract the list of datasets in the curated remote repository and focus on those concerning ENCODE:

dataset_list <- show_datasets_list(remote_url)
list <- unlist(lapply(dataset_list[["datasets"]], function(x) x$name))
grep(pattern = 'ENCODE', x = list, value = TRUE)

Select ChIP-seq data from the ENCODE NARROW dataset AUG_2017 and ENCODE BROAD dataset AUG_2017, aligned to HG19:

Enc_Broad <- read_gmql("public.HG19_ENCODE_BROAD_AUG_2017",
                          is_local = FALSE)

Enc_Narrow <- read_gmql("public.HG19_ENCODE_NARROW_AUG_2017",
                           is_local = FALSE)

Select Encode ChIP-seq aligned to HG19 and related to human embryonic stem cell line (H1-hESC):

HM_TF_rep_broad <- filter(Enc_Broad, assay == "ChIP-seq" & file_status == 
                          "released" & biosample_term_name == "H1-hESC" &                            
                          output_type == "peaks")
HM_TF_rep_narrow <- filter(Enc_Narrow, assay == "ChIP-seq" & file_status ==
                          "released" & biosample_term_name == "H1-hESC" & 
                          output_type == "optimal idr thresholded peaks")

Take the union of the two previously generated datasets:

HM_TF_rep <- union(HM_TF_rep_broad, HM_TF_rep_narrow)

Filter out samples subjected to pharmacological treatment or with specific "audit" marker:

HM_TF_rep_good_0 <- filter(HM_TF_rep, !biosample_treatments == "*" & ! (audit_error == "extremely low read depth" | audit_error == "extremely low read length") & !(audit_warning == "insufficient read depth") & !(audit_not_compliant == "insufficient read depth" | audit_not_compliant =="insufficient replicate concordance" | audit_not_compliant == "missing input control" | audit_not_compliant == "severe bottlenecking" | audit_not_compliant == "unreplicated experiment"))

Filter out samples related to HM histone modifications:

TF_rep_good_0 <- filter(HM_TF_rep_good_0, !(experiment_target == "H2AFZhuman" | experiment_target == "H3F3A-human" | experiment_target == "H3K27ac-human" | experiment_target == "H3K27me3-human" |        experiment_target == "H3K36me3-human" | experiment_target == "H3K4me1-human" | experiment_target ==   "H3K4me2-human" | experiment_target == "H3K4me3-human" | experiment_target == "H3K79me2-human" |      experiment_target == "H3K9ac-human" |  experiment_target == "H3K9me1-human" | experiment_target == "H3K9me2-human" | experiment_target == "H3K9me3-human" | experiment_target == "H4K20me1-human"))

Update region attributes with length of each region and, for each sample, compute the number of regions and the sum of each region length just created:

TF_rep_good_1 <- select(TF_rep_good_0, regions_update = list(length = right - left))

TF_rep_good <- extend(TF_rep_good_1, Region_number = COUNT(),
                     sum_length = SUM("length"))

(1) Run the processing to obtain a dataset on which finding the filtering threshold;

TF_rep_good_merged <- aggregate(TF_rep_good, groupBy =
                     conds(default = c("biosample_term_name")))


TF_rep_good_ordered <- arrange(TF_rep_good_merged,                
                      regions_ordering = list(ASC("length")))

collect(TF_rep_good_ordered, name = "TF_rep_good_ordered")
job <- execute()

(1.1)Monitor the job status:

trace_job(remote_url , job$id)

(2) Run the processing to obtain a dataset on which finding the filtering threshold; once the job status is 'SUCCESS' download the resulting dataset obtained remotely in the working directory of the local File System:

dataset_name <- job$datasets[[1]]$name
print(dataset_name)

GRL_TF_rep_good_ordered <- download_as_GRangesList(remote_url, dataset_name)

download_dataset(remote_url, datasetName = dataset_name, path = './Results_use_case_3')

# Uncomment to load previously obtained results:
#setwd('./Results_use_case_3')
#dataset_name <- '_20220105_160442_TF_rep_good_ordered' #specify the name of the already downloaded dataset
#GRL_TF_rep_good_ordered <- import_gmql(dataset_name, is_gtf = TRUE)

(3) Processing to obtain a dataset on which finding the filtering threshold; once the needed dataset is saved, perform in the R environment all the remaining operations required to retrieve the threshold value:

name_sample <- names(GRL_TF_rep_good_ordered)
g <- GRL_TF_rep_good_ordered[[name_sample]]
Region_number_tot <- length(g)  
n_up <- Region_number_tot * 0.95 
n_up_1 <- n_up + 1
index <- which(g$order >= ceiling(n_up) & g$order <= floor(n_up_1))
region <- g[index]

(4) Processing to obtain a dataset on which finding the filtering threshold; the chosen threshold is the length of the 95percentile of the region number:

threshold <- region$length
threshold <- as.numeric(threshold)
threshold

Going back to RGQML remote processing, take only the regions with region lengths greater than 100 and smaller than the threshold:

TF_rep_good_filtered_0 <- filter(TF_rep_good, r_predicate = length >= 100 & length <= threshold)

Create new metadata for each sample, with number of filtered regions and the sum of their lengths:

TF_rep_good_filtered <- extend(TF_rep_good_filtered_0,
                            region_number_filtered = COUNT(),  
                            sum_length_filtered = SUM("length"))

Combine multiple replicate samples of the same TF experiment:

TF_0 <- cover(TF_rep_good_filtered, 1, ANY(), groupBy = 
               conds("experiment_target"))

Add new region attribute as length of each region after sample combination:

TF_1 <- select(TF_0, regions_update = list(length_cov = right - left))

Create new metadata for each sample, i.e. number of combined regions and min, max and sum of their lengths:

TF <- extend(TF_1, region_number_cover = COUNT(), sum_length_cover =
              SUM("length_cov"), min_length_cover = MIN("length_cov"), max_length_cover = MAX("length_cov"))

Materialize TF dataset into repository and download it on mass memory but also in main memory as GRangesList

collect(TF, name= "TF_res")
res <- execute()

#Monitor job status:
trace_job(remote_url, res$id)

res_name <- res$datasets[[1]]$name
download_dataset(remote_url, res_name, path = './Results_use_case_3')

samples <- download_as_GRangesList(remote_url, res_name)

# Uncomment to load previously obtained results:
#setwd('./Results_use_case_3')
#dataset <- '_20220105_180300_TF_res' #specify the name of the already downloaded dataset
#samples <- import_gmql(dataset, is_gtf = TRUE)

Log out from remote engine:

logout_gmql(remote_url)

Post-processing before the analysis with TFHAZ

TF=vector()
len_rep <- sapply(samples, function(x) len <- length(x))

TF_rep <- mapply(function(x, l){
  exp <- x$experiment_target
  TF_temp <- rep(exp, l)
  TF_0 <- append(TF, TF_temp)
}, samples@metadata, len_rep) 

TF <- unlist(TF_rep)
H1_hESC_0 <- unlist(samples)
data1 <- data.frame(H1_hESC_0)
data <- cbind(data1[1:5], TF)
GR_H1_hESC <- as(data, "GRanges")

After loading the TFHAZ package, find the transcription factor HOT DNA zones focusing on one chromosome at at time (e.g. chr21), by executing the following instructions:

library(TFHAZ)
TF_acc_21_w_0 <- accumulation(GR_H1_hESC, "TF", "chr21", 0)
plot_accumulation(TF_acc_21_w_0)
d_zones <- high_accumulation_zones(TF_acc_21_w_0, method =
                                     "overlaps", threshold = "std", plotZones = TRUE)  
print(d_zones)